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ABSTRACT 

The origin and structure of the magnetic fields in the interstellar medium of spiral galaxies 
is investigated with 3D, non-ideal, compressible MHD simulations, including stratification in 
the galactic gravity field, differential rotation and radiative cooling. A rectangular domain, 
1x1x2 kpc 3 in size, spans both sides of the galactic mid-plane. Supernova explosions drive 
transonic turbulence. A seed magnetic field grows exponentially to reach a statistically steady 
state within 1.6 Gyr. Following Germano (1992) we use volume averaging with a Gaussian 
kernel to separate magnetic field into a mean field and fluctuations. Such averaging does not 
satisfy all Reynolds rules, yet allows a formulation of mean-field theory. The mean field thus 
obtained varies in both space and time. Growth rates differ for the mean-field and fluctuating 
field and there is clear scale separation between the two elements, whose integral scales are 
about 0.7 kpc and 0.3 kpc, respectively. 

Key words: galaxies: ISM — ISM: kinematics and dynamics — turbulence — MHD — 
dynamo 



1 INTRODUCTION 

The interstellar medium (ISM) of spiral galaxies is strongly af- 
fected by energy injection from supernovae (SNe), which drive 
highly compressible, transonic turbulent motions. This makes it ex- 
tremely inhomogeneous, yet it supports magnetic fields at a global 
scale of a few kiloparsecs. Mean-field dynamo models have proved 
successful in modelling galactic magnetic fields and offer a use- 
ful framework to study them and to interpret their observations 
(e.g., Beck et al. 1996; Shukurov 2007). Turbulent dynamo action 
involves two distinct mechanisms. The fluctuation dynamo relies 
solely on the random nature of the fluid flow to produce random 
magnetic fields at scales smaller than the integral scale of the ran- 
dom motions. The mean-field dynamo produces magnetic field at 
a scale significantly larger than the integral scale, and requires ro- 
tation and stratification to do so. For any dynamo mechanism, it is 
important to distinguish the kinematic stage when magnetic field 
grows exponentially as it is too weak to affect fluid motions, and 
the nonlinear stage when the growth is saturated, and the system 
settles to a statistically steady state. 

The scale of the mean field produced by the dynamo is con- 
trolled by the properties of the fluid flow. For example, in the 
simplest a 2 -dynamo in a homogeneous, infinite domain, the most 
rapidly growing mode of the mean magnetic field has scale of order 
47r?7t/a, where a can be understood as the helical part of the ran- 
dom velocity and rjt is the turbulent magnetic diffusivity (Sokoloff 
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et al. 1983). For any given a and rjt, this scale is finite and the mean 
field produced is, of course, not uniform. 

In galaxies, the helical turbulent motions and differential ro- 
tation drive the so-called aco -dynamo, where the mean field has 
a radial scale of order Ar ~ 3|D|~ 1/3 (/lR) 1/2 at the kinematic 
stage (Starchenko & Shukurov 1989), with V the dynamo num- 
ber, h c± 0.5 kpc the half-thickness of the dynamo-active layer and 
R c± 3 kpc the scale of the radial variation of the local dynamo 
number. For \V\ — 20, this yields Ar ~ 1.3 kpc. 

These estimates refer to the most rapidly growing mode of the 
mean magnetic field in the kinematic dynamo; it can be accompa- 
nied by higher modes that have a more complicated structure. Mag- 
netic field in the saturated state can be even more inhomogeneous 
due to the local nature of dynamo saturation. The mean magnetic 
field can have a nontrivial, three-dimensional spatial structure, and 
any analysis of global magnetic structures must start with the sepa- 
ration of the mean and random (fluctuating) parts. However, many 
numerical studies of mean-field dynamos define the mean magnetic 
field as a uniform field obtained by averaging over the whole vol- 
ume available, or in the case of fields showing non-trivial varia- 
tions in a certain direction, as planar averages, e.g., over horizontal 
planes for systems that show z-dependent fields (e.g., Brandenburg 
& Subramanian 2005). 

The mean and random magnetic fields are assumed to be sep- 
arated by a scale, A, of order the integral scale of the random mo- 
tions, Zo; A is not necessarily precisely equal to lo, however, and 
must be determined separately for specific dynamo systems. The 
leading large-scale dynamo eigenmodes themselves have extended 
Fourier spectra, both at the kinematic stage and after distortions by 
the dynamo nonlinearity. Thus, both the mean and random mag- 
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netic fields are expected to have a broad range of scales, and their 
spectra can overlap in wavenumber space. Thus, it is important to 
develop a procedure to isolate a mean magnetic field without un- 
physical constraints on its spectral content. This problem is espe- 
cially demanding in the multi-phase ISM, where the extreme inho- 
mogeneity of the system can complicate the spatial structure of the 
mean magnetic field. 

The definition of the mean field as a horizontal average may 
be appropriate in simplified numerical models where the vertical 
component of the mean magnetic field, (B z ), vanishes because 
of periodic boundary conditions applied in x and y; otherwise, 
V • (B) = cannot be ensured. An alternative averaging procedure 
that retains three-dimensional spatial structure within the averaged 
quantities is volume averaging with a kernel Gg (r — r'), where £ is 
the averaging length: (f) £ = f f(r')Ge(r — r f ) d 3 r f , for a ran- 
dom field /. A difficulty with volume averaging, appreciated early 
in the development of turbulence theory, is that it does not obey the 
Reynolds rules of the mean (unless £ — > oo), ((f) i g) i / (f)^(g)^ 
and ((f) e ) e / (f) e (Sect. 3.1 in Monin & Yaglom 2007). Horizon- 
tal averaging represents a special case with £ —> oo in two dimen- 
sions, and thus satisfies the Reynolds rules; however, the associated 
loss of a large part of the spatial structure of the mean field lim- 
its its usefulness. Germano (1992) suggested a consistent approach 
to volume averaging which does not rely on the Reynolds rules. 
A clear, systematic discussion of these ideas is provided by Eyink 
(2012, Chapter 2), and an example of their application can be found 
in Eyink (2005). The averaged Navier-Stokes and induction equa- 
tions remain unaltered, independent of the averaging used, if the 
mean Reynolds stresses and the mean electromotive force are de- 
fined in an appropriate, generalized way. The equations for the fluc- 
tuations naturally change, and care must be taken for their correct 
formulation. An important advantage of averaging with a Gaussian 
kernel (Gaussian smoothing) is its similarity to astronomical ob- 
servations, where such smoothing arises from the finite width of a 
Gaussian beam, or is applied during data reduction. 

Here we analyze magnetic field B produced by the rotational 
shear and random motions in the numerical model of the SN-driven 
ISM presented by Gent et al. (2012, hereafter, Paper I), using Gaus- 
sian smoothing. We suggest an approach to determine the appro- 
priate length £, and then obtain the mean magnetic field Bi. The 
procedure ensures that V • Bg — 0. The random magnetic field b 
is then obtained as b = B — Bt>. The Fourier spectra of the mean 
and random magnetic fields overlap in wavenumber space, but their 
maxima are well separated. 



2 BASIC EQUATIONS 

We model the ISM, with parameters generally typical of the Solar 
neighbourhood, using a three-dimensional Cartesian grid in a re- 
gion 1.024 x 1.024 x 2.176 kpc 3 in size, with 1.024 kpc in the 
radial and azimuthal directions and 1.088 kpc vertically on either 
side of the galactic plane. D — 1.024 kpc is the largest effective 
length scale resolved in all directions. The injection of thermal and 
kinetic energy by SNe heats the gas and drives intense random mo- 
tions, whose correlation length is of order Z — 100 pc within 
±200 pc of the galactic plane (Paper I). This part of the compu- 
tational domain encompasses of order 400 turbulent cells, so statis- 
tical properties of the ISM can be determined reliably. 

Our code (based on the Pencil Code, http://code.google. 
com/p/pencil-code/) and model have been tested to ensure that the 
relevant physical processes, including the expansion of individual 



SN remnants, are modelled reliably (Paper I). The minimum nu- 
merical grid spacing required for that is A = 4 pc, and we use 
256 x 256 x 544 mesh points (excluding boundary zones), with 
the x, y, z coordinates corresponding to the r, 0, z coordinates of 
a rotating cylindrical reference frame whose z-axis is aligned with 
the angular velocity of galactic rotation and whose x-axis points 
away from the galactic centre. 

The basic equations include the mass conservation equation, 
the Navier-Stokes equation, the heat equation — all presented in 
Paper I, but now including additional terms due to the interaction 
with the magnetic field — and the induction equation written in 
terms of the vector potential A (in the gauge $ = 77 V • A): 

^ = (U + u) x B + (n + ( V )V 2 A + (V • A)V( V , (1) 

where u is the deviation of the gas velocity from the background 
shear profile, B is the magnetic flux density, and 77 is the magnetic 
diffusivity. The advective derivative, D/Dt = d/dt+(U + u)-V, 
includes a large-scale flow U = (0, Sx, 0) due to galactic differen- 
tial rotation with the shear rate S — rdQ/dr. The angular velocity 
in the Solar vicinity is = — S = 25 km s _1 kpc -1 , but here we 
consider a model with Q = 2Q = —S to enhance the mean-field 
dynamo action (Gressel et al. 2008) and thereby save computation 
time. The system is driven by localized injection of kinetic and ther- 
mal energy modelling SN explosions. To control widespread strong 
shocks, we include shock-capturing transport coefficients, includ- 
ing (r], which differ from zero in converging flows. We use here a 
power-law parameterization of radiative cooling labelled WSW in 
Table 1 of Paper I. 

We apply periodic boundary conditions in the azimuthal (y) 
direction. Differential rotation is modelled using the shearing- sheet 
approximation with sliding periodic boundary conditions in x. 
Boundary conditions on the top and bottom faces allow gas to es- 
cape without preventing inward flows, and for the magnetic field 
we adopt dA x /dz dA y /dz = A z = ^> B x = B y = 0. 
Vertical magnetic energy flux vanishes and periodic boundaries re- 
strict the planar averages of B z to be zero. The volume or planar 
averages of B x and B y are unrestricted and may generate uniform 
contributions to the horizontal magnetic field. 

The initial conditions, described in Paper I, are close to hydro- 
static and thermal equilibrium. We include a weak initial azimuthal 
magnetic field with B y = 0.03 /iG at the mid-plane, decreasing 
with z in proportion to gas density p, so that the initial field av- 
eraged over the whole domain is about 0.1 nG. In 200-400 Myr, 
the hydrodynamic parameters of the system settle into a quasi- 
stationary turbulent state independent on the initial conditions. 



3 THE MEAN MAGNETIC FIELD 

Figure 1 demonstrates the nearly exponential growth of mag- 
netic energy density against the (approximately) stationary back- 
ground of turbulent motions and thermal structure. As the aver- 
age magnetic energy density (cb) = (\B\ 2 )/8tt becomes com- 
parable to the average kinetic energy density (e u ) at t > 1 Gyr, 
the latter shows a modest reduction, as expected for the conver- 
sion of the kinetic energy of random motions to magnetic en- 
ergy. For t > 1.4 Gyr, magnetic field settles to a statistically 
steady state (es) ~ 2.5 x 10 -13 erg cm -3 , somewhat larger than 
(e u ) ^ 1.6 x 10~ 13 erg cm -3 , while the thermal energy density 
(ex) ~ 10~ 12 erg cm -3 is not affected by magnetic field. 

Magnetic field B can be decomposed into Bz, the part aver- 
aged over the length scale £, and the complementary fluctuations 
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Figure 1. Evolution of energy densities averaged over the whole computa- 
tional domain; thermal (red), kinetic (e u ) (blue) and magnetic (c-q) 
(black). Guide line (purple, dotted) indicates exp(7.5£ Gyr -1 ). 
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Figure 2. (a) Energy densities at t = 1.2 Gyr of mean, (e^) (black, 
solid), and fluctuating, (e&) (blue, dashed), magnetic fields as functions of 
the averaging length £, normalized by ( (ee)^) and averaged over the region 
\z\ < 0.5 kpc; also (e Bg ) att = 0.8 Gyr (dotted) and t = 1.6 Gyr (dash- 
three-dotted). Values of £ sampled are indicated by red crosses, (b) Deriva- 
tive of (e&) from Panel (a) with respect to £, with the same line types. 

b, 

B = B e + b, B e = (B) v b = B-B e , (2) 
using volume averaging with a Gaussian kernel: 

(B) e (x)= f B(x')G t (x-x')d 3 x', (3) 
Jv 

Ge(x) = (27r^ 2 )" 3/2 exp [-x 2 /{2£ 2 )} , 

where V is the volume of the computational domain. This opera- 
tion preserves the solenoidality of both Bi and b and retains their 
three-dimensional structure. For computational efficiency, the aver- 
aging was performed in the Fourier space where the convolution of 
Eq. (3) reduces to the product of Fourier transforms. 

Since the averaging (3) does not obey the Reynolds rules for 



the mean, the definitions of various averaged quantities should be 
generalized as suggested by Germano (1992). In particular, the lo- 
cal energy density of the fluctuation field is given by 

e b {x) = ^-[ \B(x')-B £ (x)\ 2 Gi(x-x')d 3 x'. (4) 

OTT Jy 

This ensures that (es) £ = es t + e&, where es = |^| 2 /(8tt) 
and eB £ = \Bg\ 2 /(8tt). Expanding B(x') in a Taylor series 
around x and using f Ge(x — sc') d 3 x' = 1 (normalization) and 
j v (x—x')Gi(x—x')d 3 x' = (symmetry of the kernel), it can be 
shown that e b = |6| 2 /(87r) + 0(£/L) 2 , where L is the scale of the 
averaged magnetic field defined as L 2 = \Bn |/| V 2 Bi | in terms of 
the characteristic magnitude of Bg and its second derivatives. Thus 
the difference between the ensemble and volume averages rapidly 
decreases as the averaging length increases, £/L — »• 0, or if the av- 
eraging is performed over the whole space, L — >> oo; this quantifies 
the deviations from the Reynolds rules for finite £/L. 

The appropriate choice of £ is not obvious. We consider a 
range < £ < 500 pc, applying the averaging to 37 snapshots 
of the magnetic field between t = 0.8 and 1.7 Gyr; the results for 
t — 0.8, 1.2 and 1.6 Gyr are shown in Fig. 2. The smaller is £, the 
closer the correspondence between the averaged field and the orig- 
inal field (since the average is effectively sampling a smaller local 
volume), and hence the smaller the part of the total field consid- 
ered as the fluctuation. Hence (B 2 ) is a monotonically decreasing 
function of £, whereas (b 2 ) monotonically increases (Fig. 2a). Our 
aim is to identify that value of £ where the variation of Bg and b 
with £ becomes weak enough. To facilitate this, we consider the rate 
of change of the relevant quantities with £, shown in Fig. 2b (note 
the different scale of the horizontal axis in this panel). The length 
£ w 50 pc is clearly distinguished: all the curves in Panel (b) are 
rather featureless for £ > 50 pc. The values in Fig. 2 have been 
obtained from the part of the domain with \z\ < 0.5 kpc, where 
most of the gas resides and where dynamo action is expected to 
be most intense. The results, however, remain quite similar if the 
whole computational domain \z\ ^ 1 kpc is used. While this value 
of £ has been estimated in a rather heuristic and approximate man- 
ner, the analysis of the magnetic power spectra in Section 4 con- 
firms that £ = 50 pc is close to the optimal choice. 

To put the estimate £ w 50 pc into context, we note that it is 
about half the integral scale of the random motions, lo, in a similar 
model of Paper I. Simulations of a single SN remnant in an ambient 
density similar to that at z = 0, reported in Paper I, also show that 
the expansion speed of the remnant reduces to the ambient speed 
of sound when its radius is 50-70 pc. 

Figure 3 illustrates the total, mean and random magnetic fields 
thus obtained. For this saturated state, the field has a very strong 
uniform azimuthal component and a weaker radial component. The 
orientation of the field is the same above and below the mid-plane 
(By > and B x < 0), with maxima located at \z\ ~ 0.2 kpc; the 
results will be reported in detail elsewhere. 



4 SCALE SEPARATION 

Scale separation between the mean and random magnetic fields 
in natural and simulated turbulent dynamos remains a controver- 
sial topic. The signature of scale separation sought for is a pro- 
nounced minimum in the magnetic power spectrum at an interme- 
diate scale, larger than the energy-range scale of the random flow 
and smaller than the size of the computational domain. The spec- 
tra of the mean and random magnetic fields obtained by Gaussian 
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Figure 3. Field lines of (a) the total magnetic field B, (b) its averaged part Bg, (c) the fluctuations b, obtained by averaging with £ = 50 pc, for t = 1.625 Gyr. 
Field directions are indicated by arrows. The colour of the field lines indicates the field strength (colour bar on the left), whereas the vectors are coloured 
according to the strength of the azimuthal (y) component (colour bar on the right). 
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Figure 4. Power spectra of total (black, dotted), mean (red, dash-dotted) 
and fluctuating (blue, dashed) magnetic fields, for Gaussian smoothing with 
£ = 50 pc, at t = 1.05 Gyr. Vertical red (dash-dotted) line marks the 
averaging wavenumber 2tt£~ 1 and vertical black (dotted) line indicates the 
Nyquist wave number ttA~ 1 . Short vertical segments mark the energy- 
range scales 2-irLJ 1 and 27rL b ~ 1 of the mean and random magnetic fields, 
respectively, obtained from Eq. (5). 



smoothing, shown in Fig. 4, have maxima at significantly different 
wavenumber s. We note that the spectrum of the total magnetic field 
does not have any noticeable local minima and, with the standard 
approach (e.g., based on horizontal averages), the system would be 
considered to lack scale separation. 

Figure 4 shows the magnetic power spectra of B, Bi and b. 
Note that £ is not located between the maxima in the power spec- 
tra of Bi and b; in fact, the spectral density of the mean field is 
negligible for k ~ 2tt£~ 1 . This can be understood from the trans- 
form of the kernel Gi(x), i.e. Gi(k) — ex.p(—£ 2 k 2 /2); this ker- 
nel would divide a purely sinusoidal field equally into the mean and 
random parts at the wavelength A eq = y/2/ In 2tt£. For £ — 50 pc, 
A eq = 0.27 kpc, and the latter figure is a better guide to the ex- 
pected separation scale. The separation of scales is immediately 
apparent in the spectra of the mean and random fields, Mi{k) and 
Mb (/c), with the former having a broad absolute maximum at about 
0.56 kpc and the latter, a broad maximum near 0.2 kpc. The ef- 
fective separation scale A « 0.48 kpc can be identified where 
Mi(k) = Mb(k), i.e. where the curves cross at log k c± 1.1. 

The integral scales of Bg and b can be obtained from their 
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Figure 5. As Fig. 4, but for (a) £ = 20 pc and (b) £ = 100 pc. 



spectral densities as 

r>7r/A 



7T r* i 

= - / /c _1 M(/c)d/c 

2 J2tt/D 



r/A 



M(k)dk 



2-k/D 



(5) 



where A is the numerical grid separation and D is the size of the 
domain (Section 12.1 Monin & Yaglom 2007). This yields Li ~ 
0.67 kpc and Lb ~ 0.28 kpc for Bg and 6, respectively (Fig. 4). 

As the magnetic field strength grows, the magnitudes of the 
spectral densities change but the characteristic wavenumbers vary 
rather weakly. Li remains close to 0.7 kpc throughout the kine- 
matic phase, increasing only beyond 1.1 Gyr to about 0.9 kpc. Lb 
increases from 0.23 to 0.28 kpc during the kinematic phase, but 
rises to 0.4 kpc after the system saturates. The stability of Li is 
consistent with the eigenmode amplification of the mean magnetic 
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Figure 6. Evolution of magnetic energy densities, averaged for full domain: 
total magnetic field (black, dotted), mean (black, solid), and random (blue, 
dashed) obtained from Gaussian smoothing with £ = 50 pc. Energy densi- 
ties of mean and random magnetic fields obtained from horizontal averaging 
are shown long-dashed (black and blue, respectively). 



Table 1. Exponential growth rates of energy in the mean and random mag- 
netic fields, T and 7, respectively, with associated values of reduced x 2 . 
From exponential fits to the corresponding curves in Fig. 6, for 0.8 < t < 
1.05 Gyr. (Growth for t < 0.8 Gyr is similar to this interval - see Fig. 1.) 
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field as expected for a kinematic dynamo, and supports t « 50 pc 
as a reasonable choice of the averaging scale. 

Figure 5 presents magnetic energy spectra obtained using 
£ = 20 pc and 100 pc. In the former, the effective separation scale 
A « 0.16 kpc is less than Lb c± 0.21 kpc. This is inconsistent, 
implying that energy at scales about Lb lies predominantly within 
Bi. For £ — 100 pc, A ^ 1 kpc is greater than Li c± 0.81 kpc, 
which is also inconsistent. Significantly, applying £ — 50 pc satis- 
fies Lb < A < Lg. Hence the scale A, at which the dominant en- 
ergy contribution switches between the mean and fluctuating parts, 
is here consistent only with £ ~ 50 pc. 

As well as different spatial scales, the mean and random mag- 
netic field energies have different exponential growth rates (Fig. 6), 
given in Table 1 . Results from horizontal averaging are also given; 
they differ from those obtained with Gaussian smoothing, espe- 
cially for the mean field. The growth rate F of the mean-field en- 
ergy is controlled by the shear rate, mean helicity of the random 
flow, and the turbulent magnetic diffusivity. Any mean magnetic 
field is accompanied by a random field, which is part of the mean- 
field dynamo mechanism (as any large-scale magnetic field is tan- 
gled by the random flow); the energy of this part of the small-scale 
field should grow at the same rate F as the mean energy. The fluc- 
tuation dynamo produces another part of the random field whose 
energy growth rate depends on the turbulent kinematic time scale 
l/u, and the magnetic Reynolds and Mach numbers. The differ- 
ence between F and 7 obtained suggests that both mean-field and 
fluctuation dynamos are present in our model. The growth rate of 
the mean magnetic energy is roughly double that of the random 
field for both types of averaging. This is opposite to what is usu- 
ally expected, plausibly because of the inhibition of the fluctuation 
dynamo by the strongly compressible nature of the flow. 



5 DISCUSSION 

The approach used here to identify the appropriate averaging length 
£ (and thus the effective separation scale A) is simple (but not over- 
simplified); £ can in fact depend on position, and it can remain 
constant in time only at the kinematic stage of all the dynamo ef- 
fects involved. Wavelet filtering may prove to be more efficient than 
Gaussian smoothing in assessing the variations of £. 

At the kinematic stage, the spectral maximum of the mean 
field is already close to the size of the computational domain, and it 
cannot be excluded that the latter is too small to accommodate the 
most rapidly growing dynamo mode. These results should there- 
fore be considered as preliminary with respect to the mean field; 
simulations in a bigger domain are clearly needed. 

The physically motivated averaging procedure used here, pro- 
ducing a mean field with three-dimensional structure, may facili- 
tate fruitful comparison of numerical simulations with theory and 
observations. Although Gaussian smoothing does not obey all the 
Reynolds rules, it is possible that a consistent mean-field theory can 
be developed, e.g. in the framework of the r-approximation (see 
e.g. Brandenburg & Subramanian 2005, and references therein). 
This approach does not rely upon solving the equations for the 
fluctuating fields, and hence only requires the linearity of the av- 
erage and its commutation with the derivatives. The properly iso- 
lated mean field is likely to exhibit different spatial and temporal 
behaviour than the lower-dimensional magnetic field obtained by 
two-dimensional averaging. 
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